Two routes to the one-dimensional discrete nonpolynomial Schrodinger equation 
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The Bose-Einstein condensate (BEC), confined in a combination of the cigar-shaped trap and 
axial optical lattice, is studied in the framework of two models described by two versions of the one- 
dimensional (ID) discrete nonpolynomial Schrodinger equation (NPSE). Both models are derived 
from the three-dimensional Gross-Pitaevskii equation (3D GPE). To produce "model 1" (which 
was derived in recent works), the 3D GPE is first reduced to the ID continual NPSE, which is 
subsequently discretized. "Model 2", that was not considered before, is derived by first discretizing 
' the 3D GPE, which is followed by the reduction of the dimension. The two models seem very 

different; in particular, model 1 is represented by a single discrete equation for the ID wave function, 
while model 2 includes an additional equation for the transverse width. Nevertheless, numerical 
analyses show similar behaviors of fundamental unstaggered solitons in both systems, as concerns 
their existence region and stability limits. Both models admit the collapse of the localized modes, 
£f) ■ reproducing the fundamental property of the self-attractive BEC confined in tight traps. Thus, 

£N| ' we conclude that the fundamental properties of discrete solitons predicted for the strongly trapped 

self-attracting BEC are reliable, as the two distinct models produce them in a nearly identical 
form. However, a difference between the models is found too, as strongly pinned (very narrow) 
discrete solitons, which were previously found in model 1, are not generated by model 2 - in fact, in 
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Mh ' agreement with the continual ID NPSE, which does not have such solutions either. In that respect, 

(-h , the newly derived model provides for a more accurate approximation for the trapped BEC. 

PACS numbers: 03.75.Lm; 05.45.Yv 

The dynamics of a dilute quantum gas which forms the Bose-Einstein condensate (BEC) is very 
accurately described by the three-dimensional Gross-Pitaevskii equation (3D GPE). This equation 
treats effects of collisions between atoms in the condensate in the mean-field approximation. In 
experimentally relevant settings, the BEC is always confined by a trapping potential. In many cases, 
the trap is designed to have the "cigar-shaped" form, allowing an effective reduction of the dimension 
• ■ from 3 to 1. In turn, the ID dynamics of the trapped condensate may be controlled by means of 
an additional periodic potential, induced by an optical lattice (OL), which acts along the axis of 
■ the "cigar". If the OL potential is sufficiently strong, the eventual dynamical model reduces to a 
ID discrete equation. In both the continual and discrete versions of the ID description, a crucially 
important feature is the form of the nonlinearity in the respective equations. In the limit of low density, 
the nonlinearity is cubic the same as in the underlying 3D GPE. In the general case, a consistent 
derivation, which starts from the cubic nonlinearity in 3D, leads to ID equations with a nonpolynomial 
nonlinearity, the respective model being called the "nonpolynomial Schrodinger equation" (NPSE). 
The discrete limit of the latter equation, corresponding to the action of the strong axial OL potential, 
was derived and investigated recently. An essential asset of both versions of the NPSE, continual 
and discrete ones, is that they predict the onset of the collapse (formation of a singularity in the 
condensate with attraction between atoms) in the framework of the ID description, thus complying 
with the fundamental property of the BEC which was predicted by the underlying 3D GPE and 
observed experimentally. However, in the case when the OL potential is very strong, an alternative 
way to derive the ID discrete model may start with the discretization of the 3D GPE, followed by the 
reduction of the dimension in the cigar-shaped trap. In this work, we report a new discrete model 
("model 2") derived in this way, which seems very different from the previously known discrete ID 
NPSE (which we call "model 1"). In particular, while model 1 amounts to a single discrete equation for 
the ID complex wave function, model 2 incorporates an additional equation for the transverse width. 
Nevertheless, numerical analysis performed in the present work shows remarkably similar behavior 
of fundamental localized modes, in the form of unstaggered discrete solitons, in both systems. The 
similarity pertains to the existence region for the solitons and their stability limits. Importantly, 
both models admit the collapse, and produce similar predictions for the collapse threshold. Thus, 
basic properties of discrete solitons found in the model of the strongly trapped BEC are trustworthy, 
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as they are reproduced independently by the two very different models. Nevertheless, a difference 
between the two models is also found: very narrow discrete solitons, which exist in model 1, are absent 
in model 2. In fact, the continual ID NPSE does not give rise to such extremely narrow solutions 
either, thus indicating that the newly derived model, although having a more complex mathematical 
form, eventually provides a more accurate approximation. 



I. INTRODUCTION 



The dynamics of Bose-Einstein condensates (BECs) made of dilute ultracold gases of bosonic atoms obeys the 3D 
Gross-Pitaevskii equation (GPE) in the mean- field approximation An effective ID equation can be derived from 
the 3D GPE to describe the BEC dynamics in prolate ("cigar-shaped") traps In the simplest case, which 

corresponds to a sufficiently low BEC density, the reduced ID equation amounts to the cubic nonlinear Schrodinger 
equation (NLSE) 0. A significant restriction in the use of the ID cubic NLSE in this context is its failure to predict 
the onset of the collapse of localized modes, which was theoretically predicted in 3D models and experimentally 
observed in the self-attractive BEC [13, [H| ■ This problem can be resolved by the more accurate reduction of the 3D 
GPE to the ID nonlinear Schrodinger equation with a nonpolynomial nonlinearity (NPSE) Q, without imposing 
the constraint of a very low density. The ID NPSE model with self-attractive nonlinearity enables the description of 
the c ollap se dynamics, and yields results which are accurately reproduced by direct simulations of the underlying 3D 
GPE [12]. An intermediate approximation, which can be obtained from the expansion of the NPSE, is represented by 
the ID NLSE with the self-focusing cubic and quintic terms. This approximation may be sufficient for some particular 
purposes [3j, [y, [7[ . 

On the other hand, the BEC trapped in a very deep optical lattice (OL) can be well described by the corresponding 
discrete equations. In particular, discrete forms of the ID GPE with the cubic nonlinearity [HI, 0, EH and ID NPSE 
[HI , have been studied in detail. Basic features of the ID continual models describing the BEC trapped in a deep 
OL find their counterparts in the discrete models, a significant one being the ability of the discrete ID NPSE to 
describe the onset of the collapse predicted by the corresponding continual equation [ijj • 

Thus, the previously explored ID discrete BEC models relied upon the discretization after the reduction of the 
dimension from 3 to 1 within the framework of the continual equations. However, an alternative approach is possible 
too, in the situation when the OL is stronger than the cigar-shaped potential: one should first discretize the 3D GPE 
in the axial direction, and then reduce the dimension to 1 @, opposite to the previously developed derivation 0]. 
The purpose of the present work is to derive the discrete model in this way, which, to the best of our knowledge, 
has not been done before. Naturally, the two different routes lead to quite different discrete systems. In particular, 
the one derived in this work involves two sets of discrete variables (the wave function and transverse width), rather 
than the single set on which previously studied models were based. Our purpose is not only to derive the alternative 
system, but also to study generic properties of fundamental unstaggered solitons in it, and compare the results with 
those obtained in the earlier known discrete NPSE. In fact, we will conclude that, despite a very different form of 
the new model, it gives rise to results quite similar to those produced by the discrete NPSE. Thus, despite significant 
differences in the form of the "competing" models, the physical predictions for the discrete BEC solitons are essentially 
the same, which attests to the reliability of these results. 

The paper is structured as follows. The new model is derived in section 2. Basic results for the existence, stability 
and dynamical behavior of fundamental discrete solitons in it are reported and compared to their counterparts 
predicted by the usual ID discrete NPSE in section 3. The paper is concluded by section 4. 



II. THE DERIVATION OF THE DISCRETE ONE-DIMENSIONAL SYSTEMS 



4] 



The starting point of the consideration is the 3D GPE, which governs the evolution of macroscopic wave function 
(r, t) of the dilute BEC near zero temperature [lj . This equation can be derived from the following action functional 
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-(x 2 +y 2 )-i 7 (A-l)H 2 UM), (1) 



where 7 = 4irT?a s /m is the strength of the interaction between bosons, a s the s-wave scattering length, N the 
number of bosonic atoms, m the atomic mass, and u>± the frequency of the transverse harmonic confinement. The 
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corresponding 3D GPE is written as 



-— V 2 + Vbcos (2kz) 
2m 
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(2) 



Expressing the length in units of a±, time in units of a;^ 1 and energy in units of Tu)±, the 3D GPE is cast in the 
scaled form, 
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~(y± + j^) + V cos (2kz) + \{x 2 + y 2 ) + 2rf>| 2 
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where T = 2(N — l)a 8 /a± is the effective strength of the self-interaction, and a± — y/h/ (muj±) is the characteristic 
length of the transverse confinement. Notice that a s < 0, i.e., T < in Eq. ([3]), corresponds to the attraction between 
bosons. 

Below, we consider the derivation of two alternative forms of the ID discrete approximation. The previously known 
"model 1" is derived starting with reduction of the 3D GPE to the ID NPSE, which is then discretized p~6l. fl7l| . Here, 
we briefly recapitulate this route of the derivation, for the purpose of the comparison with novel "model 2" , which is 
obtained from the 3D GPE by first discretizing it, and then reducing the resulting system to the ID form. 



A. Model 1: the reduction of the dimension followed by the discretization 



The derivation of model 1 starts with the dimensional reduction, which is performed by the minimization of action 
functional ([1]) , choosing the wave function in the form of 



„ exp[-(x 2 +;/ 2 )/(2a 2 (z,t))] 



(4) 



where real a is the local width of the transverse confinement, and complex ID (axial) wave function / is normalized 



\m\'dz = \ (5) 

Inserting trial wave function Q in Eq. (TTJ) and performing the integration over x and y, the action can be written 



as 



S = 



(0) 



in the approximation which neglects da/dz Q. The Euler-Lagrange equations derived from Eq. © by varying with 
respect to /* and a, 
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v 4 = i + r|/| 2 , 

may be combined into the equation (NPSE) derived in Ref. namely, 
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Assuming that OL potential in this equation is strong enough, one can further derive the discrete version of the ID 
NPSE [IB) 113] • To this end, the continual wave function is approximated by a superposition of orthonormal modes 
W n (such as Wannier functions), which are tightly confined in a vicinity of local potential minima, z n = nn/q, with 
integer n: 



f(z,t)=J2fn(t)W n (z), 
n 



(10) 
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/„ being the respective complex amplitudes. This decomposition is made unique by imposing a condition that the 
largest value of each local function W n is 1 [HI]. Next, one may insert this ansatz in Eq. ([7]), multiply the resulting 
equation by the complex conjugate of the local mode, and integrate over z. The so derived ID discrete NPSE, together 
with the discrete version of equation (JSJ), take the form of 
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where the local norms and parameters e, C and g are 

1 d 2 
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In the tight-binding approximation, C is positive definite. 

In what follows below, the discrete system based on Eqs. (fTTj) and (TT2^) is referred to as "model 1" . These equations 
correspond to the Lagrangian, 
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/n + C/n(/n+l + /n-l) 



Model 1 conserves two dynamical invariants, viz., the norm (alias "power") and Hamiltonian, 

= £ |/„| 2 , H = £ ~ /™+l| 2 + Vl-9\fn\ 2 \fnf 



(14) 
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B. Model 2: the discretization followed by the reduction of the dimension 



The derivation of the novel system starts with the direct discretization of the 3D GPE, i.e., Eq. @, by adopting 
the following ansatz for the wave function, 



ip(x,y,z,t) = ^2<l> n (x,y)W n (z) ) 



(16) 



where W n (z) is the same set of local modes as in Eq. (|T0|) . We insert this ansatz in Eq. ([2]), multiply the resulting 
equation by W*, and integrate the result over z. This leads to the semi-discrete equation, 
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where acts on coordinates x and y, and parameters e, C and g are defined as per Eqs. {T 
Equation (|17p can be derived from the corresponding Lagrangian, 
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The further simplification is performed by substituting ansatz 

1 



n (x,y,t) 



• exp 



x 2 + y 2 



fn(t) 



(17) 



(18) 



(19) 



into Lagrangian (|18[) and integrating over the (x, y) plane, which eventually leads to an effective Lagrangian for ID 
discrete fields f n (t) and a n (t), 

d 1/1 



L = J2 Wn 
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The Euler-Lagrange equations derived from Eq. ([20)) by varying the Lagrangian with respect to /* and a n are 
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dt In ~ 
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0. 



Equations (I21[) and (|22f conserve the respective norm and Hamiltonian, which are [cf. Eq. I|15p] 

^ = £i/„i 2 , 



(22) 
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+ fn-lfn) + 7T~2"l.fn| 



(23) 



It is worthy to note that Eq. (|12p in model 1 involves only er rl (without coupling to cr„±i), because the corresponding 
spatial derivative, da/dz, was neglected in the underlying expression ([6]). On the contrary to that, Eq. ([22]) in model 2 
couples a n to a n ±i, as this coupling is implicitly retained by expressions (fT5|) and (JTSJ) . For that reason, cr 2 cannot be 
explicitly found from Eq. (f2"2"|) , hence model 2, unlike model 1, cannot be reduced to a single equation for the discrete 
wave function. Nevertheless, model 2 is amenable to a consistent numerical analysis, see below. The two models can 
be made formally equivalent only if transverse widths a n are postulated to be constant (t- and n-independent), in 
which case either system reduces to the standard discrete ID NLSE. Finally, it is worthy to mention that both models, 
i.e., Eqs. HJ), (fT2"l) . on the one side and Eq. (|2"T]) on the other, take identical limit forms in the anticontinuum limit, 
C = 0. 



III. FUNDAMENTAL BRIGHT SOLITONS 
A. The existence of fundamental solitons 



In this section we present families of fundamental bright-soliton solutions to Eqs. (f2Tj) and ([22]) . identify their 
stability and compare them to the corresponding soliton families in model 1 [l6l [l7| . If the nonlinearity is attractive 
(g < 0), the stationary solutions of models 1 and 2 with chemical potential are sought for by the substitution into 
Eqs. (TTJ), |H2J) and fl2TJ), respectively, of 



fn(t) = \g\ 1/2 «n exp {-i/ri) . 



(24) 



In doing so, we assume that e = 2C was fixed by means of an obvious additional shift of the chemical potential. Then, 
real discrete functions u n obey stationary equations, which are 



in model 1, and 
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in model 2. In the latter case, widths a n are to be found from Eq. ([22]) . which takes the form of 
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FIG. 1: Values of the transverse width, a n , as functions of the chemical potential, /i, for fundamental on-site solitons in model 1 
(the ordinary model) are shown by curves which are marked by triangles for the central site, full circles - for the first neighbors, 
and empty circles - for the second neighbors to the central site. In model 2 (the new system) the values of a„ for fundamental 
on-site solitons are shown by dashed lines for the central site, dotted lines - the first neighbors, and dashed-dotted lines - the 
second neighbors to the central site. The inter-site coupling constant is C — 0.2 (a) and C = 0.8 (b). 




FIG. 2: The same as in Fig. [T] but for families of inter-site solitons found in both models 1 and 2. 



Stationary equations (l25l) and (|26[) were solved numerically, using an algorithm based on the modified Powell 
minimization method [14 , Il6j ] . The initial ansatz used to construct on-site and inter-site-centered discrete solitons in 

model 1 was, respectively, = (..., 0, A, 0, ...) and (..., 0, A, A, 0, ...), where A is a real constant obtained from 

Eq. |25|) in the corresponding approximation. These soliton solutions are then used as an initial ansatz to generate 
discrete solitons in model 2. Results reported below were obtained in the lattice composed of 101 or 100 sites, for the 
on-site and inter-site configurations, respectively. It was checked that the results do not alter if a larger lattice had 
been used. 

We start the presentation of the results by plotting, in Figs. [1] and transverse widths a n versus the chemical 
potential for stationary unstaggered solitons of the on-site and inter-site types, respectively, found in both models 
1 and 2 [recall that the vanishing field corresponds to a n — 1, see Eqs. ([T2|) and (j22|) ]. It is evident that these 
characteristics of the soliton families are close for both models, with some difference observed in central parts of the 
solitons. 

To outline the entire existence region for the fundamental solitons, we followed the usual approach, identifying it as 
the region where CW (continuous-wave) solutions are modulationally unstable, hence the existence of bright solitons 
should be expected. The CW solutions can be easily found from stationary equations ([23]) and (|2"rj)) in the form of 

u n = Ue-^, ficw = [1 + (3/2) gU 2 ] /^(1 + gU 2 ), a 2 n = 1 + gU 2 , (28) 

which shows that in the case of the attractive interaction, g = —1, the amplitude of the CW solution is subject to 
constraint U < 1, and the respective chemical potential takes values /j,qw < 1- Then, straightforward calculations 
yield a dispersion relation for frequency Q and wavenumber q of small modulational perturbations around the CW 
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FIG. 3: The CW solutions are unstable to modulational perturbations in regions below curves fJ.{q), which are determined by 
condition Q 2 > 0, see Eq. (|29[) . The solid and dashed curves correspond to models 1 and 2, respectively. The figure pertains 
to (7 = 0.8. 



solution: 



n 2 = -2A 



2A 



2gU 2 



9 2 U A 



V 1 + f}U 2 2(1 + gU 2 ) ^l+gU 2 + D) 



(29) 



where A = 2Csin 2 (q/2), and D = A in model 2, D = in model 1. Analysis of this dispersion relation demonstrates 
that, in either model, the CW solution is unstable for /x < 1 against long- wave perturbations (for small values of q), 
see Fig. ©. 

In the case of the repulsive contact interactions, g > 0, one may expect the existence of staggered discrete solitons 
[lH . The analysis of the modulational instability of the respective staggered CW solutions (not given here in detail) 
shows that, for both models, the staggered CW states are indeed unstable in their entire existence region (/x > 1), 
which indicates that staggered solitons may exist at /x > 1 in either model. The further study of model 1 indicates 
that it supports stronglylocalized (tightly pinned) staggered solitons of the fundamental type (as reported before, 
see Figs. 20-23 in Ref. 16]). Model 2 reproduces only low-amplitude tightly pinned staggered modes near the lower 
boundary of the existence region. Those solutions are not displayed here, as the objective of the work is to focus on 
unstaggered solitons in the case of the self-attraction, when the difference of both models, 1 and 2, from the ordinary 
discrete cubic NLSE is most essential. 



B. The norm and free energy of fundamental unstaggered solitons 

The power (norm) of unstaggered solitons in the two models, P — Yln=i l w n| 2 > an( i their free energy G = H — /xP, 
where H is obtained by inserting stationary solution into expressions (|2"2l or (f2"B"|). are compared in Figs. ([!]), 
(0 and ^ for two values of the intersite-coupling constant, C — 0.2 and C = 0.8. These values, which were used 
in Figs. [1] and OS actually correspond, severally, to limit cases of strongly discrete and quasi-continual systems. It is 
seen that the norms and free energies of the soliton families in both models feature similar behaviors: the same sign 
of the slope, dP/d^i, for both the on-site and inter-site modes, the equality of the on-site and inter-site free energies 
in the region of a small power and small |/x| (wide solitons with a small amplitude). Quantitative differences between 
the norm and free energy between the two models are larger for higher values of C, i.e., closer to the continuum limit. 



8 




FIG. 4: The norm ("power") P versus chemical potential /x for on-site fundamental solitons, in both models 1 and 2. The 
coupling constants are C = 0.2 (a) and C = 0.8 (b), which correspond to strongly discrete systems and quasi-continual systems, 
respectively. 
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FIG. 5: The norm ("power") P versus chemical potential /i for inter-site fundamental solitons, in models 1 and 2. The coupling 
constants are C = 0.2 (a) and C = 0.8 (b). 



Some differences also appear in the region of large norms, which corresponds to very narrow solitons, where the P(fi) 
curves show a trend to saturation. In particular, in that region amplitude of the narrow solitons in models 2 may 
exceed 1, which is the absolute upper limit (collapse threshold) in model 1, as follows from Eq. ([8]). In fact, numerical 
results indicate at the existence of a similar limit in model 2, although it does not explicitly follow from Eqs. (|2~lj) 
and (|22[1 . Finding exact values for this limit is complicated by difficulties in obtaining accurate numerical results for 
very narrow (strongly pinned) solitons in model 2. 




FIG. 6: Free energy G versus norm P for on-site and inter-site fundamental solitons, in both models 1 and 2, for C = 0.2 
(a) and C = 0.8 (b). Solid and dashed lines correspond to the on-site solitons in models 2 and 1, respectively. Dotted and 
dashed-dotted lines represent, respectively, inter-site solitons in models 2 and 1. 
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FIG. 7: Pure real unstable eigenvalues ("ev") for inter-site unstaggered solitons in models 1 and 2 are shown by black and 
white circles, respectively, for C = 0.2. Note that, in this case, pure real eigenvalues for on-site solitons have not been found. 



The sufficient condition for soliton stability in the NPSE is the spectral condition according to which the corre- 
sponding eigenvalues, found from linearized equations for small perturbations, must not have a positive real part 
[T3. [l9|. Although the Vakhitov-Kolokolov (alias slope) stability criterion, dP/dfi < is not strictly applicable 
to the models with nonpolynomial nonlinearity we will also briefly comment it. 

In Figs. (j4|) and ((5]), P(/i) curves for the fundamental-soliton families have the same sign of the slope (dP/dfi) in 
models 1 and 2. Therefore, in the region of small P and small positive fi, the on-site solitons may be stable in both 
models, according to the VK criterion. In all other cases, the on-site solitons are expected to be unstable. On the 
other hand, the VK criterion indicates the stability of inter-site solitons in the entire existence region, in both models. 

To confirm the stability of the solitons, we have checked the spectral condition within the framework of the linear 
stability analysis, following the approach elaborated in Refs. [HI, [ljj])- The results show that the fundamental 
unstaggered solitons in both models indeed have the same stability properties, as seen in Figs. ([7]) and ©. In 
particular, for on-site solitons the stability takes place only in a narrow region in the plane of ([i, C). This region is 
characterized by very close values of the norms of the on-site and inter-site solitons, hence a very small difference in 
the respective values the free energy, which, in turn, implies a very low Peierls-Nabarro barrier (lfjl. [Tt|. i.e., mobility 
of the solitons. 

All conclusions concerning the mobility, perturbed evolution, and development of the collapse instability are quali- 
tatively identical in both models. In particular, unstable inter-site modes which find on-site counterparts with close 
values of the norm (or free energy), and unstable on-site modes with close inter-site counterparts evolve into robust 
breathers, with almost no loss of the norm. In the opposite case, the instability leads to collapse of the localized 
modes. Our calculations show that collapse thresholds are nearly identical in models 1 and 2. In Ref. [16| . it was 
found, in the framework of model 1, that (as mentioned above) in soliton solutions cannot exceed 1; actually, 
— 1 in Eq. (|25|) is the collapse threshold, which is determined by the singular structure of the on-site nonlinearity. 
A similar trend is demonstrated by numerical results obtained in model 2. 

Summarizing this section, we conclude that both discrete models, 1 and 2, translate fundamental properties of 
bright solitons, known in the ID continual NPSE, into the discrete setting in essentially the same way (fo r the discrete 
model of type 1, the correspondence with the continual equation was demonstrated in recent work [2fj). The most 
significant among these properties is the possibility to describe the onset of the collapse in the framework of the ID 
geometry. 



The mean-field dynamics of BEC confined in a combination of a cigar-shaped trap and axial OL (optical-lattice) 
potential, is well approximated by the ID NPSE (nonpolynomial Schrodinger equation), in both continual [J, |20( and 
discrete settings [is, 0|. In particular, both versions of the NPSE admit a possibility to describe the onset of the 
collapse of the localized modes in the framework of the ID approximation, thus making this approximation compliant 
with the well-known fundamental property of the self-attractive BEC confined in tight traps. The objective of the 
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FIG. 8: The same as in Fig. [7] but at C = 0.8. In this case, plots (a) and (b) display the eigenvalues for the on-site solitons 
and inter-site solitons, respectively. 

present work was to verify the relevance of the description based on the discrete limit of the NPSE-type equation, by 
comparing two alternative versions of this approximation: the previously derived one ("model 1"), which is based on 
the reduction of the underlying 3D GPE to ID continual NPSE, followed by its discretization, and the new "model 
2" , whose derivation starts with the discretization of the 3D GPE, followed by the reduction of the dimension. The 
two steps of the derivation, viz., the dimension reduction and discretization, are apparently non-commutative, and, 
accordingly, final forms of the two models appear to be vastly different; in particular, model 1 reduces to a single 
discrete equation for the wave function, while model 2 includes an additional equation for the transverse width. 
Nevertheless, the numerical solutions for families of fundamental unstaggered bight solitons produce very similar 
results in both models, as concerns their existence and stability, and, especially, the crucially important feature - the 
collapse threshold. A difference between the two models was found too, as model 2 (the newly introduced one) does 
not reproduce the region of high-amplitude strongly pinned narrow unstaggered and staggered modes, which was 
predicted in model 1. Actually, that feature of model 1 is impugnable, as it is not reproduced by the continual NPSE 
with the strong OL potential, as noted in Ref. [T^. In this respect, model 2 (the novel one), although having a more 
cumbersome mathematical form, seems physically preferable. 
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